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AN EFFICIENT MESHFREE IMPLICIT FILTER FOR NONLINEAR 

FILTERING PROBLEMS* 

FENG BAQt, YANZHAO CAO*, CLAYTON G. WEBSTERt, AND GUANNAN ZHANGt 

Abstract. In this paper, we propose a meshfree approximation method for the implicit filter 
developed in which is a novel numerical algorithm for nonlinear filtering problems. The implicit 
filter approximates conditional distributions in the optimal filter over a deterministic state space grid 
and is developed from samples of the current state obtained by solving the state equation implicitly. 
The purpose of the meshfree approximation is to improve the efficiency of the implicit filter in 
moderately high-dimensional problems. The construction of the algorithm includes generation of 
random state space points and a meshfree interpolation method. Numerical experiments show the 
effectiveness and efficiency of our algorithm. 
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1. Introduction. Nonlinear filters are important tools for dynamical data as¬ 
similation with applications in a variety of research areas, including biology [T, 20 


mathematical finance J4 


multi-target tracking 15 


11 , signal processing 14 20,23 , image processing 22 , and 


18 . To put it succinctly, nonlinear hltering is an extension 


of the Bayesian framework to the estimation and prediction of nonlinear stochastic 
dynamics. In this effort, we consider the following nonlinear hltering model 


^ = fit, Xt, Wt), (state) 

Yt = g{t,Xt) + Vt, (observation) 


( 1 . 1 ) 


where / and g are two nonlinear functions, {Xt G > 0} and {Yt G > 0} 
are the stochastic state and observation processes, respectively, {Wt G W,t > 0} is 
a random vector representing the uncertainty in Xt, and {Vt G > 0} denotes 

the random measurement error in Yt. In the discrete setting, the nonlinear hltering 


model in (1.1) takes the form 


Xk = fk(Xk-i,Wk-i), 
Yk = gk(Xk) + Vk, 


(state) 

(observation) 


( 1 . 2 ) 


where G M'" and G K® are mutually independent white noises. 

Let Yi:fc := cr{Fi,T 2 ,- • • ,Lfc} denote the cr hied generated by the observational data 
up to the step k. The goal of nonlinear hltering is to hnd the posterior probability 
density function (PDF) of the state Xk, given the observation data Yi.,k, so as to 
compute the quantity of interest (Qol), given by 

E[$(Afc)|Ti,fc] = inf {E[|$(Afc) - Zj^]; Z G Zk} , 
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where $(•) is a test function, and Zk denotes the space of all Yi:fe-measurable and 
square integrable random variables. 


Tremendous efforts have been made to solve nonlinear filtering problems in the 
last few decades. Two of the well-known Bayesian filters are extended Kalman filters 
(EKFs) [S pOpbpTp^ , and particle filters 5|[7pp3 . The key ingredient of the EKFs 
is the linearization of both / and g in (1.11, so that the standard Kalman filter can 
be applied directly. Thus, if the nonlinearity of the state and the observation systems 
is not severe, then the EKFs can provide efficient and reasonable inferences about 
the state, otherwise, the performance of the EKFs can be very poor. For particle 
filters, the central theme is to approximate the desired posterior PDF of the state 
by the empirical distribution of a set of adaptively selected random samples (referred 
to as “particles”). The particle filter method is essentially a sequential Monte Carlo 
approach, which requires no assumption on the linearity of the underlying system. 
As such, with sufficiently large number of particles, it is capable of providing an 
accurate approximation of the posterior PDF for a highly nonlinear filtering problems. 
However, there are some fundamental issues concerning the efficiency and robustness 
of particle filters [^. For example, since the empirical PDF is constructed based on 
particles with equal weights after resampling, the particle filter still needs a lot of 
samples in order to accurately approximate the target distribution. 


To overcome such a disadvantage, the authors proposed a new nonlinear filter 
named “implicit filter” . This approach adopts the framework of Bayesian fil¬ 
tering, which has two stages at each time step, i.e., prediction and update. At the 
prediction stage, we estimate the prior PDF p{Xk\Yi.k-i) of the future state Xk given 
the current available observation information Yi-k-i] at the update stage, we update 
the prior PDF by assimilating the newly received data to obtain the estimate of 
the posterior PDF p{Xk\Yi.k). The implicit filter is distinguished from the particle 
filters by the use of interpolatory approximations to the prior and posterior PDFs. 
Specifically, in the particle filter, p(Xfc|Yi:fc_i) is approximated by explicitly propagat¬ 
ing the samples of the current state Xk-i\Yi:k-i through the nonlinear state equation 
Yk = /fc(Alfc_i, Wfc_i), and constructing the empirical PDF of Xk\Yi.k-i. In the im¬ 
plicit filter, the interpolation of p{Xk\Yi:k-i) requires its function values at a set of 
grid points of the future state Xk- Under the condition that Xk = a; S we solve 
implicitly the state equation x = fk{Xk-i,Wk-i) given a set of Monte Carlo samples 
of Wk-i, so that the value of p{Xk = x\Yi.k-i), at the grid point of x, can be esti¬ 
mated by averaging the function values of p{Xk-x\Yi.k-i) at all the solutions of the 
state equation. As an implicit scheme, the implicit filter has a stabilizing effect which 
provides more accurate numerical approximations to the solution of the nonlinear 
filtering problem than the particle filter method [^. 


The main challenge of the implicit filter method is that the conditional PDF of 
the nonlinear filtering solution is estimated at grid points. As such the method suffers 
the so called “the curse of dimensionality” when the dimension of the state variable is 
high. In addition, the efficiency of the method may be significantly reduced when the 
domain of the PDF is unbounded. In this paper, we propose to construct a meshfree 
implicit filter algorithm to alleviate the aforementioned challenges. Motivated by the 
particle filter method, we first generate a set of random particles and propagate these 
particles through the system model and use these particles to replace the grid points in 
the state space. After that we generate other necessary points through the Shepard’s 
method which constructs the interpolant by the weighted average of the values on 
state points 12 . In order to prevent particle degeneracy in the generation of random 











A meshfree implicit filter for nonlinear filtering problems 


3 


state points, we introduce a resample step in the particle propagation. In addition 
we choose state points according to the system state, which make them adaptively 
located in the high probability region of the PDF of state. In this way, we solve 
the nonlinear filtering problem in a relatively small region in the state space at each 
time step and approximate the solution on a set of meshfree state points distributed 
adaptively to the desired PDF of the state. Furthermore, since we approximate the 
PDF as a function on each state point, instead of using state points themselves to 
describe the empirical distribution, the implicit filter algorithm requires much fewer 
points than the particle filter method to depict the PDF of the state. 

The rest of this paper is organized as follows. In we introduce the mathemat¬ 
ical framework of the Bayesian optimal filter. In ^ we construct meshfree implicit 
algorithm. In ^ we demonstrate the efficiency and accuracy of our algorithm through 
numerical experiments. Finally, contains conclusions and directions for the future 
research. 

2. Bayesian optimal filter. For m,n G N’*', let Xm-.n and Fmin denote the a 
fields generated by {Xm, X^+i, ■ ■ •, Xn} and {Ym, Ym+i, ■■■, Y^}, respectively. For 
k = N+, we use Xk to represent a realization of the random variable Xk, and define 


p{xk\-) ■■= p{Xk = Xk\-) 


for notational simplicity. It is easy to see that the dynamical model in (1.2) is Marko¬ 
vian in the sense that 


p(^Xk\Xi-k—ljYi-k — l') — p{xk\Xk—l')- 

We also know that the measurements are conditionally independent given Xk, i.e., 

p{Yk\Xi-_kjYi-_k—l) — p(Yk\xk^’ 

The Bayesian optimal filter constructs the conditional distribution p{xk\Yi.k) recur¬ 
sively in two stages: prediction stage and update stage. 

For k = 1,2,---, assume that p{xk-i\Yi:k-i) is given. In the prediction stage 
p{xk\Yi:k-i) is evaluated through the Chapman-Kolmogorov formula: 


p(a:fe|yi:fe_i) = / p{xk\xk-i)p{xk-i\Yi,k-i)dxk-i. (2.1) 


In the update stage, the prior PDF obtained in (2.1) is used to obtain the posterior 
PDF p(xfc|Yi:fc) via the Bayes’ formula: 


p{xk\Yi,k) = 


piYk\xk)p{xk\Yi,k-i) 

piYk\Yi,k-i) 


p(Yk\xk)p{xk\Yi,k-i) 

\xk)pixk |Fl:fc—l) dXk 


( 2 . 2 ) 


3. The meshfree implicit filter. In this section, we construct the meshfree 
implicit filter algorithm. The algorithm is based the implicit filter algorithm on grid 
points [^. The implicit filter algorithm introduced in is developed from the gen¬ 
eral framework of the Bayesian optimal filter discussed above, in which the primary 
computational challenge is the numerical approximation of the term p(xk\xk-i) in 


( 2 . 1 ). 
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3.1. The prediction stage. For k = 1,2, - , the goal of this stage is to ap¬ 

proximate the prior distribution p{xk\Yi.}~_i) of the state X}., given the posterior 
distribution p(a:fe_i|yi:fc-i) of the state Xk-i- Due to the the fact that 


p{Xk\Xk-l) = ^w[p{^k\Xk-l->'Wk-l)] = / p{Xk\Xk-l,Wk-l) ■ p{Wk-l)dWk-l, 
the prior PDF p(xfc|Yi:fc_i) derived in identity (2.1) can be rewritten as 

p{xk\Yi,k-i) = Ei^[p{xk\xk-i,Wk-i)]p{xk-i\Yi,k-i)dxk-i, (3.1) 

JR‘‘ 


where £„,[•] represents the expectation with respect to the white noise Wk-i, and the 
PDF p{xk\xk-i,Wk-i) is 


p{Xk\Xk-l,Wk-l) 


OO^ Xk — fk{xk — l: Wk—l') ^ 
O 5 ^ fk{xk — ljWk—l\ 


(3.2) 


with p{xk\xk-i,Wk-i)dxk = 1 for any Xk-i G and Wk-i G The def¬ 
inition in ( |3.2[ ) can be viewed as a generalization of the Dirac delta function in 
the space R^x R'^ x R’', where the mass is located according to the state equation 
Xk — fk{xk—l') Wk—l') • 

Note that the estimation of (3.1) requires an approximation to the expectation 
Eyj[p{xk\xk-i,Wk-i)\- To this end, we first draw M independent samples 
of the white noise Wk-i, and define an approximation to p{xk\xk-i,Wk-i) as 


M 

Tr^{xk\xk-i,Wk-i) (xk\xk-i,Wk-i), (3.3) 

i=i 


with 


S^j^^^{xk\xk-i,Wk-i) 


00 , Wk-i = wl_^ and Xk = fk{xk-i,wl_^), 
0, otherwise, 


which is essentially a restriction oip{xk\xk-i,Wk-i) in the subset There¬ 

fore, the expectation E^[p{xk\xk-i,Wk-i)] in (3.1) can be approximated by 


E»[p(a:fe|a;fc-i, Wfc-i)] « Eu, {xk\xk-i,Wk-i)] , 

M „ 

= Y/ {xk\xk-i,Wk-i)p{wk-i)dwk-i. 

^ Jr- ''-1 


(3.4) 


To construct an interpolation of p{xk\Yi:k-i), the next step is to approximate 
p(a;fc|yi:fc-i) at a point set Tik ■= {x\}fLi C R"^ with N G N"*". By substituting 


Xk = x]. into (3.1)-(3.4), we have 


p{.xl\Yi,k-i) = [ Eyj[TT^{xl\xk-i,Wk-i)]p{xk-i\Yi.,k-i)dxk-i+'Jllik-i, (3-5) 
Jr‘^ 


where7?.^|j,_;^ := p{xl\Yi,k-i)-{xl\xk-i,Wk-i)]p{xk-i\Yi,k-i)dxk-i is the 
approximation error. Then, by further fixing Wk-i = w^k-i-! location of the mass 
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of S j in the space of Xk-i, denoted by can be obtained by 

implicitly solving the state equation 

h (4-i> 


w 


k-1 


= xl, j = !,■■■ ,M, 


which is the reason we named the approach the implicit hlter. Now substituting x\. 
into (3.4), and using the same sample set as above to approximate the 

integral on the right hand side of (3.4), we obtain 

M 


E, 




[iT^ {xl\xk-l,Wk-l)] = 5] I ^ 


M 




M 


(3.6) 


= — y (5 3 
M ^ K-i 






then replacing ¥.-uj\7t^ {xWxk-i^Wk-i)] in (3.5) with (3.6), we have 

M 


pixl\Yi,k-i) = 


/R'i 

1 

M 


1 

M 


i=i 


p{xk-i\Yi,k-i)dxk-i + 








(3.7) 

where p(a;^’;ij^|Fi:fe_i) is the value ofp (a;fe_i|li:fc_i) dX x''^^_^. Neglecting the error term 

(3.7), we obtain the following iterative numerical scheme for constructing 
an approximation, denoted by p(a;^|yi:fc_i), of the prior PDF p(a;yFi:fe_i), i.e., 


g{x\\Yi.,k-i) = ^ ^ Q{xll^\Yi,k-i). (3.8) 

i=i 

In our previous work [^, the subsets "Hs,, for fc = 0,1,..., were defined by a full 
tensor product mesh, denoted by 

M := X X ■ ■ ■ (3.9) 

on a d-dimensional hyper-cube [oi, 6i] x • • • x [ad, 6d], where = 1,..., d, is 

a uniform partition of the interval with grid points. It is simple to 

implement but has several significant disadvantages. First, at each time step, one 
needs to approximate the prior PDF p(a;fe|yi:fe_i) at a total of x • • • x N^d) 
points which grows exponentially as the dimension d increases. This is also known 
as “the curse of dimensionality”. On the other hand, since the construction of M. 
is not informed by the target PDF, the domain [ai,5i] x • • • x [ad,hd] needs to be 
defined sufficiently large, so as to capture the statistically significant region of the 
PDF. This may lead to a great waste of computation effort in the low probability 
region of p(a;fc|Ti:fc_i). 

To alleviate such disadvantages, we propose to develop a distribution-informed 
meshfree interpolation approach to efficiently approximate the prior PDF. The central 
idea of the generation of random points for the state variable is to build a set of points, 
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denoted by Hk, according to the state distribution. To begin with, we generate 
Hq = of N random samples from the initial PDF pq of the initial state: 

■^0 := with xf, = f. 


If the initial PDF pq is close to the true state distribution, it’s obvious that our 
random state points are more concentrated near the target state. For k = 1,2, ■ ■ ■ , K, 


we propagate points to through the state equation (1.2|: 




.iV, 


where are N random samples according to the PDF of Wk-i- Denote 

Hk ■= and approximate the conditional PDF pixklYi-.k-i) on Hk with the 

scheme given by (3.81. In this way, the random points in Hk move according to the 


state model. As opposed to particle filter methods, which use the number of particles 
to represent empirical distributions and require a large number of particles to follow 
the state distribution, in the implicit filter method we provide an approximation of 
the value of the PDF at each state point. Therefore, much fewer points are needed to 
describe the state PDF and the random state points are not necessary to accurately 
follow the state distribution. 

3.2. The update stage. By incorporating the new data Yk, we update the prior 
PDF p{xk\Yi:k-i) at each grid point x\, using the Bayesian formula, to obtain 


p{xl,\Yi,k) = ^piYk\xl)p{xl\Yi,k-i) 

G/c 

= ^PiYk\xl)g{xl\Yi,k-i)+TZl^k> 


(3.10) 


where Q{xWYi,k-i) is given in (3.8), Ck is the normalization factor, and := 

d^p(Xk\xi){p{xl\Yi,k-i) - g{xl\Yi,k-i)) is the approximation error. By neglecting 
the error term in (3.10), we obtain the following iterative numerical scheme for 


the update stage on Hk, i-e.. 


9{xl\Yi:k) = —p{Yk\xl)g{xl:\Yi,k-i) 


(3.11) 


where g{xl.\Yi,k) is desired the approximation of the posterior PDF p{x].\Yi:k). 

Next, we use interpolation methods to construct the approximation g{xk\Yi,k) of 
p{xk\Yi,k) from values {g{xWYi,k)}xi^Hk 


gixk\yi:k) = 9i^k\yi-.k)(l>’‘{xk), 


(3.12) 


xieT-ik 


where is the set of basis functions. Since the state points in Hk are generated 

randomly in the meshfree framework, standard polynomial interpolation is unstable 
due to the uncontrollable Lebesgue constant. Instead, we propose to use the Shepard’s 
method [12| , which is an efficient meshfree interpolation technique, to construct the 
interpolant g{xk\Yi,k)- The basic idea of the Shepard’s method is to use the weighted 
average of {g{Xk\Yi,k)}xieHk interpolating approximation. Specihcally, for a 
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given point Xk & we re-order the points in Tik by the distances to Xk to get a 
sequence such that 

\\xk - < ||a;fc - if h < h, 

where || • || is the Euclidean norm in Then, for a pre-chosen integer L we use the 
first L values in {gix^k'^ \ Yi.k)}fLi to approximate g{xk\Yi.k) as follows 

L 

g{xk\Yi-.k) = E Qi^k^\yi:k) * hi{xk), (3.13) 

1^1 


where the weight hiixk) is defined by 




hi{xk):= . 

E/=i lla^fc -4 II 


Note that X]i=i hi{xk) = 1- From (3.13), we have 


Q{Xk\yi-.k) -p{Xk\yi-.k) = X! 


1^1 


+ ^p{x^k\yi-.k) ■ hi{xk) -p{xk\yi:k), 


1^1 


where 

L L 

'^pix^k'^Xi-k) ■ hi{xk) -p{xk\Yi,k) = ^ (p{xll:^\Yi:k) - pixk\Yi:k)'^ ■ hi{xk) (3.14) 
1=1 1=1 


is the error of the Shepard’s interpolation. We assume that p{xk\Yi.,k) has bounded 
first order derivative. For each pair p(a:fe| Yi:fe) andp(a;^^^ 1^1,^) the approximation error 


\p{xl!'^\Yi.,k) — p{xk\Yi-,k)\ is controlled by the distance 

(i) 




and the derivative 


p'(z\Yi-,k), where z is a point between Xk and x\ . It is reasonable to assume that 
in high probability region of the derivative p'{z\Yi.,k) is large. It’s worth pointing 
out that the random state points generated in this algorithm are concentrated in the 
high probability region. Thus, if Xk lies in the high probability region, the distance 
\\xk — x^^W is small, which balances the error brought by the large derivative. On the 
other hand, if Xk lies in the low probability region, although the distance \\xk — 
is relatively large, the approximation error (3.14) is still small due to the small value 
of the derivative p'{z\Yi.k). 


3.3. Resampling. Similar to the particle filter method, the above random state 
points generation suffers from the degeneracy problem for long term simulations, 
especially for high-dimensional problems. After several time steps, the probability 
density tends to concentrate on a few points which dramatically reduces the number 
of effective sample points in Hk- 

In this work, we propose an occasional resampling procedure to address these 
problems and rejuvenate the random points cloud. At the time step fc — 1, the re¬ 
sampling procedure takes place after we obtain g{xk-i\Yi.k-i), in order to remove 
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the degenerated points in T-Lk-i using the information provided by Q(xk-i\Yi.,k-i). 
Specifically, the first step is to develop a degeneracy metric to determine the neces¬ 
sity of doing resampling. To this end, we define the following degenerated subset 

^k—l ^ 1; 

Sk-1 = e 'Hfe-i,p(4_i|ri:fe_i) < e}, (3.15) 

where e > 0 is a user-defined threshold. We also define 




to be the index set of Sk-i- Then, the degeneracy of T-Lk-i can be measured by the 
ratio #(5fe_i)/#('Hfc_i) G [0,1], where #(•) denotes the number of points in a set. If 
the ratio is smaller than a threshold r € [0,1], then we will skip the resampling step 
and propagate "Hfc-i to get Hfe; otherwise, the set ?{k-i is considered degenerated, 
and the resampling procedure is needed. 

In resampling, instead of propagating ?{k-i to "Hfe, we aim at constructing an 

intermediate point set, denoted by Hk-i ■= {s;! i }fli and propagate T-Lk-i- through 
1—1 ^ ^2 2 
the state model (1.2) to obtain Tik- According to the definition of Sk-i in (3.15), we 

consider the state points in Tik-i\Sk-i are in the statistically significant region of 
Q(xk-i\Yi.k-i), so that we first put those points in i.e.. 


for iiJ{Sk-i). 

1% 2 

For the state points in Sk-i, we replace them by generating new samples from 
Q{xk-i\Yi:k-i) using the importance sampling [^, i.e., 

xl_i ~ gixk-i\Yi,k-i) for i G J{Sk-i). 

As a result, the resampling procedure helps us remove the state points with low 
probabilities, and makes the state point set Tlk concentrated in the high probability 
region of the posterior PDF p(xfe_i |Yi:fc_i) at each time step. 


3.4. Summary of the algorithm. Finally, we summarize the entire meshfree 
implicit filter algorithm introduced in §3.1^j|T3|in Algorithm 1 below. 


4. Numerical experiments. In this section, we present two numerical exam¬ 
ples to examine the performance of our meshfree implicit filter method. In Example 1, 
we use a two dimensional nonlinear filtering problem to show the distributions of the 
random points Tik- In Example 2, we solve a three dimensional bearing-only tracking 
problem, which is a six dimensional nonlinear filtering problem. For this higher di¬ 
mensional problem, we compare the accuracy and efficiency of our meshfree implicit 
filter method with the extended Kalman filter and the particle filter. 


Example 1. In this example, we consider the two dimensional noise perturbed 


tumoral growth model 21 


dXt = F{Xt)dt + a ■ dWt, (4.1) 

where Wt is a two dimensional standard Brownian motion and cr = (0.01,0.01)^. The 
state process Xt = (A(^,Af)’^ is a two dimensional vector, F{Xt) := {fi{Xt), f 2 {Xt))'^ 
is defined as 
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Algorithm 1: The meshfree implicit filter algorithm 


Initialization: set the number of samples M for estimating the number 

of state points N, the resampling thresholds e and r 
while fc = 1, 2, • • • , do 

Compute the ratio 4f{Sk-i)/ff{TLk-i) 
if #(5fc_i)/#(H/c-i) < T then 


Propagate Hk-i through the state model (1.2 1 to obtain TLk 
else 

Resample and construct the intermediate state set TLk-i 


Propagate through the state model ( |1.2[ ) to obtain Tik 

end if 


Prediction: solve Q{xk\Yi.,k-i) using (3.8|, at each point in TLk 


Update: solve Q{xk\Yi.,k) using (3.11) and (3.13) 

end while 


and 


h{Xt) = a^X} - asXf ■ (X^)- 


Here, fi models the Gompertzian growth rate of the tumor and /2 gives the degree 
of vascularization of the tumor which is also called “ angiogenic capacity”. 

To approximate the state variables, we discretize the dynamic system (4.1) in 
time and obtain a discrete state model 


Xk — F{Xk-i) • A + cr • u)k-i- 


(4.2) 


Here, cok is a two dimensional zero mean Gaussian white noise process with covariance 
Q = /A, where I is the 2x2 identity matrix and A is the time partition stepsize. 
The measurement of the state model is given by 

Yk = {Xl,Xlf + R-Vk, 


where Vk is a two dimensional zero mean Gaussian white noise process with covariance 
A = I A, / is a 2 X 2 identity matrix and R = (0.1,0.!)^. 

In the numerical experiment, we use uniform time partition with stepsize A = 0.2 
and simulate the state process for if = 40 with initial state Xq = (0.8,0.3)^ and 
parameters ai = 1, a 2 — 0 . 2 , 0:3 = 0 . 2 . At time step A: = 0, we initialize the prior 
PDF po by iV(Ao)E), where Aq = (0.78,0.32)^ and 


f 0.052 Q \ 

0 0.l2 ) ■ 


(4.3) 


In Figure [4T| we plot 1500 random samples generated from the initial PDF pq, 
which are our initial random points TLq- Figure |4^ illustrates the behavior of random 
state points TLk at time steps k = 1,2,3,10,20,40, respectively. In Figure |T^ the 
blue dots in each figure plot the random state points obtained by using the dynamic 
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Fig. 4.1: Example 1 : Initial random state space points "Ho 


state points generation method introduced in Sectionj^and the red cross in each figure 
gives the true state Xk at the corresponding time step. From the figures we can see 
that all the points are moving according to the state model and are concentrated 
around the true state. 


To present the accuracy of the algorithm, we show the simulation of the tumoral 

The black curves are the true and X^ coordinate 


growth states in Figure 4.3 


values of the tumoral growth states, respectively. The blue curves show the simulated 
states obtained by using the meshfree implicit filter method. 


Example 2. In this example, we study a six dimensional target tracking problem. 
In Figure |4^ the target, denoted by the red line, moves in the three dimensional space 
and two platforms on the ground, denoted by pentagons, take angular observations 
of the moving target. 

The state process X^ = {Xl, X^, X^, X^, X^, X^)"^ is described by the following 
dynamic model 


— fiXk-i) + cr ■ uJk-i, 


(4.4) 


where {X^,X‘^, X^) describes the position of the moving target which is controlled by 
parameters {X'^,X^,X^). The system noise Wfc = (w^, is a zero 

mean Gaussian white noise process with covariance Q = /A, / is the 6x6 identity 
matrix and A is a given time period, cr = (0.1, 0.1,0.1, 0.01, 0.01,0.01)^ is a constant 
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X' 

(a) k = 1 


0.8 

0.7 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

X' 

(b) k = 2 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 

x' 

(c) k = 3 


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

x' 

(d) k = 10 




0.1 - 

0 - 
0.1 


0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

X' 

(e) k = 20 


0.1 - 
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0.1 
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(f) k = 40 


Fig. 4.2: Example 1: State space points Hk at time steps k = 1,2,3,10,20,40. 


vector and / is given by 


f{Xk) = 


+sin(aA:f_i)A 

xl_^ + 

A|_i+1^2A 

y A^_i+t;3A j 
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Fig. 4.3: Example 1: Simulation of the tumoral growth states 


The measurements Y/c of the state process from the two locations are given by 
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where Vfc is a 4 dimensional zero mean Gaussian white noise process with covariance 
A = I A, / is a 4 X 4 identity matrix, R = (0.6,0.6, 0.6,0.6)"^, (ui, 6 i) and ( 02 , 62 ) are 
locations of two observers. We choose A = 0.3, a = 3, = U 2 = U 3 = 0.05 . Also, 

we assume that platforms are located at (ai, 6 i) = (16,6), ( 02 , 62 ) = (8,15) and the 
initial sate is given by a Gaussian N{Xo, E) where = (2, 2,1,0.4, 0.4,0)^ and 

/l^OOO 0 0\ 

0 l2 0 0 0 0 

00 1^0 0 0 

^ “ 0 0 0 0.2^ 0 0 
0 0 0 0 0.2^ 0 

\ 0 0 0 0 0 0.22 y 

The target will be observed over the time period 0 < t < 15. In the numerical 
experiments, we compare the performance of our meshfree implicit filter with the 
extended Kalman filter and the particle filter. In particular, we compare the estimated 
mean values of the states process along each dimension in Figure In the particle 
filter method, we choose 15,000 particles. In the meshfree implicit filter method, 
we choose the number of state points to be TV = 4,000 and the number of random 
samples in the implicit filter Monte Garlo simulation to be M = 6. The black curves 
in Figure show the real states process along each direction, the green curves give 
the estimated means obtained by the extended Kalman filter method, the red curves 
give the estimated means obtained by the particle filter method, and the blue curves 
give the estimated means obtained by the meshfree implicit filter. We also plot the 
error err^ corresponding to all three methods in figure [4(6| As we can see from figure 


extended Kalman filter and the implicit filter is the most accurate approximation in 
this experiment. 

To further compare the efficiency between the meshfree implicit filter and the 
particle filter, we repeat the above experiment over 50 realizations and show the 
average GPU time and the corresponding global root mean square error errc defined 

by 


4.5 and|4.6| the implicit filter and the particle filter are much more accurate than the 


err^ = 


, , 50 A 

j = l k=l 


4.1 


we can see 


where errk{j) is the error of the j-realization at time step k. In table 
that with 15,000 particles, the CPU time of the particle filter method is comparable 
to that of the implicit filter with 4,000 random state points, but the global RMSE of 
the particle filter is more than doubled the RMSE of the implicit filter. With 25, 000 
particles, the particle filter method achieves an accuracy comparable to the implicit 
filter, but at a significantly higher cost. 


5. Conclusions. In this work, we proposed an efficient meshfree implicit filter 
algorithm by evaluating the conditional PDF on meshfree points in the state space. 
These meshfree points are chosen adaptively according to the system state evolution. 
We also apply Shepard’s method as the meshfree interpolation method to compute 
interplants with random state points. In order to address the degeneracy of the 
random points, we use importance sampling method to construct a resample step. 
Numerical examples demonstrate the effectiveness and efficiency of our algorithm. In 






14 


F. Bao, Y. Cao, C. G. Webster and G. Zhang 




(a) (b) 




(c) (d) 




(e) 


(f) 


Fig. 4.5: Example 2 : Comparison of estimated states, (a) Shows the comparison on 
direction, (b) Shows the comparison on X^ direction, (c) Shows the comparison 
on X^ direction, (d) Shows the comparison on X'^ direction, (e) Shows the comparison 
on X^ direction, (f) Shows the comparison on X^ direction. 
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Fig. 4.6: Example 2 : Comparison of L? error. 


Table 4.1: Example 2: Efficiency comparison 


Methods 

CPU time (seconds) 

errc 

Implicit filter (4, 000 state points ) 

83.14 

0.0924 

Particle filter (15,000 particles) 

82.89 

0.2545 

Particle filter (20, 000 particles) 

142.61 

0.1687 

Particle filter (25, 000 particles) 

209.27 

0.1057 


the future, we plan to perform a rigorous numerical analysis for the meshfree implicit 
filter algorithm. 
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